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We evaluated the C*6 coefficients of Yb-Yb and Yb -alkali/group II van der Waals interactions 
with 2% uncertainty. The only existing results for such quantities are for the Yb-Yb dimer. Our 
value, C§ — 1929(39) a.u., is in excellent agreement with the recent experimental determination of 
1932(35) a.u. We have also developed a new approach for the calculation of the dynamic correction 
to the blackbody radiation shift. We have calculated this quantity for the Yb 6s 2 1 So — 6s6p s Po 
clock transition with 3.5% uncertainty. This reduces the fractional uncertainty due to the blackbody 



radiation shift in the Yb optical clock at 300 K to 10 



level. 



PACS numbers: 06. 30. Ft, 34.20.Cf, 32.10.Dk, 31.15.ac 



I. INTRODUCTION 

Ytterbium (Yb: Z=70) has recently emerged as a sub- 
ject of great interest in ultracold chemistry, physics, and 
metrology. For example, the first state-resolved observa- 
tion of ultracold chemical reactions was recently reported 
for the Yb+ + Rb — > Yb + Rb+ system Q. Yb is a 
favorite candidate for the studies of ultracold gas mix- 
tures. For example, Li and Yb mixtures have recently 
been brought to simultaneous quantum degeneracy [2j- 
5] . Controlled production of ultracold YbRb* molecules 
by photoassociation in a mixure of Rb and Yb gases was 
recently reported in Q . Such mixtures are of interest for 
producing ultracold polar molecules for study of dipo- 
lar quantum matter, fundamental symmetry studies, and 
many-body quantum simulation [5|. The availability of 
5 bosonic and 2 fermionic stable isotopes makes Yb es- 
pecially attractive for studies of multicomponent super- 
fluids. 

The spectrum of Yb contains a number of long-lived 
excited states that arc conveniently accessed by opti- 
cal techniques. This makes Yb an excellent candidate 
for atomic parity violation (APV) studies that test the 
Standard Model of electroweak interactions, put limits 
on its possible extensions, constrain parameters of weak 
hadronic interactions, and may yield information on neu- 
tron distributions within nuclei 0,1]. The APV signal 
recently observed in the Yb 6s 2 1 S , - 5d6s 3 L>i 408- 
nm forbidden transition [8|, [9( is two orders of magni- 
tude larger than in Cs, subject of the most accurate APV 
study to date. Such long-lived states are also convenient 
for the development of next-generation ultra-precise fre- 
quency standards. The Yb ^So — > 3 Pq 578-nm transition 
now provides one of the world's most accurate optical 
atomic frequency standards [13, [H| ■ 

The work carried out in this Letter is pertinent to all 
applications mentioned above. Our two main subjects 
are the determination of van der Waals Cg coefficients 



that characterize the long-range interactions between two 
atoms, and the blackbody radiation (BBR) shifts of the 
two states in atomic clock transitions. Knowledge of the 
long range interactions in Yb-Yb and Yb- alkali/group II 
dimers is critical to understanding the physics of dilute 
gas mixtures. The dynamic correction to the BBR shift 
is one of the largest irreducible contributions to uncer- 
tainty budget of the Yb clock [10( , and it is difficult to de- 
termine experimentally. These two seemingly disparate 
topics both require accurate determination of frequency- 
dependent atomic polarizabilitics over a wide range of 
frequencies. Therefore, it is natural to consider them in 
the same work. A future accurate theoretical determi- 
nation of the APV amplitude in Yb requires a similar 
approach, and this work provides a background for such 
studies. 

We carry out the calculation of frequency-dependent 
atomic polarizabilities using the first-principles approach 
that combines configuration interaction (CI) with the 
coupled-cluster all-order approach (CI+all-order) that 
treats both core and valence correlation to all orders. 
Several new method developments are presented in this 
work. First, we have implemented the reduced linear 
equation (RLE) and direct inversion in iterative sub- 
space (DIIS) stabilizer procedures described in Ref. [l2[ 
into the coupled-cluster part of the CI+all-order method. 
Otherwise, the construction of the effective Hamiltonian 
needed for the incorporation of the core and core- valence 
correlations into the CI method could not be carried out 
due to convergence problems associated with extremely 
large correlations involving the 4f shell. Second, we have 
applied the CI+all-order method for the first time to the 
calculation of Cq coefficients. Finally, we have developed 
a new approach to the calculation of the dynamic cor- 
rection to the blackbody radiation (BBR) shift in terms 
of the second partial derivative with respect to frequency 
of the dynamic polarizability, as obtained from the solu- 
tion of the inhomogeneous equation in the valence sec- 
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tor. Previous calculations of the dynamic correction to 
the BBR shift accounted for the contributions of just a 
few intermediate states to the polarizability 13]. 

Unless stated otherwise, we use atomic units (a.u.) for 
all matrix elements and polarizabilities throughout this 
paper: the numerical values of the elementary charge, 
e, the reduced Planck constant, Ti = h/2n, and the 
electron mass, m e , are set equal to 1. The atomic 
unit for polarizability can be converted to SI units via 
a/h [Hz/(V/m) 2 ]=2.48832xl(r 8 a (a.u.), where the con- 
version coefficient is 4ireoa,Q/h, &o is the Bohr radius and 
eo is the electric constant. 

Calculation of Yb properties requires an accurate 
treatment of both core- valence and valence- valence corre- 
lations. This can be accomplished within the framework 
of the CI+all-order method that combines configuration 
interaction and coupled-cluster approaches 6j j . Here 
we report the extension of this method that resolves the 
convergence problems associated with particularly large 
correlation corrections as well as apply it for the first 
time to the calculation of the C§ coefficients. We refer 
the reader to Refs. [l4], EH] for detailed description of 
this approach, and here we report only new method de- 
velopments specific to this work. In order to establish the 
accuracy of our approach, we also perform the pure CI 
and the CI combined with many-body perturbation the- 
ory (CI+MBPT) calculations carried out with the same 
parameters such as basis set, configuration space, number 
of partial waves, etc.. 

The single-electron energies and the wave functions are 
found from the solution of the Dirac-Hartree-Fock (DHF) 
equations. Then the wave functions and the low-lying 
energy levels are determined by solving the multiparti- 
cle relativistic equation for two valence electrons fl7j . 
H c ff(E n )$ n — E n Q n . The effective Hamiltonian is de- 
fined as H c g(E) — f/pc + where Hpc is the Hamil- 
tonian in the frozen-core approximation. The energy- 
dependent operator £(-£?) which takes into account vir- 
tual core excitations is constructed using second-order 
perturbation theory in the CI+MBPT method [l8| and 
using a linearized coupled-cluster single-double method 
in the CI+all-order approac h 1141 . However, the CI+all- 
order approach developed in [l4|,ll5j could not be directly 
implemented for Yb owing to convergence problems of 
the all-order equations associated with large oscillations 
of the iterative solution due to very large correlations in 
the 4/ shell. Both the Yb 2+ core and some of the Yb + va- 
lence shell all-order equations that are used to construct 
the effective Hamiltonian diverge using conventional it- 
eration schemes. We have resolved this problem by u sing 
RLE and DIIS convergence stabilizers described in [l2[ 
within the framework of the CI+all-order method. The 
main idea of these approaches is to accumulate several 
iterations and determine a next best solution based on 
all stored data. Convergence was achieved for the ns, 
np, and (n — l)d valence states with n = 6 — 9. 

We present the energy levels obtained in the CI, 
CI+MBPT, and CI+all-order approximations and com- 



pare them with the experimental values [19J in Table I 
of the supplementary material (2(j. At the CI stage, the 
theoretical energy levels differ rather significantly from 
the experimental energies, up to 19% for the 6s6p states. 
Including the core- valence correlations in the second or- 
der of the MBPT improves the agreement to the 1.5-5.5% 
level. Further improvement of the theoretical energies is 
achieved when the CI+all-order approximation is used. 
The two-electron binding energy is accurate to 0.7% with 
the ab initio CI+all-order approach, a factor of 2 im- 
provement in comparison with the CI+MBPT result. 

The valence part of the polarizability is determined 
by solving the inhomogeneous equation of perturbation 
theory in the valence space, which is approximated as 

{E v - JlrfOI* («, M')> = D eS , q \* (v, J, M)> (1) 

for a valence state v with the total angular momentum 
J and projection M [5l[. The effective dipole operator 
D e g includes random phase approximation (RPA) cor- 
rections. The ionic core part of the polarizability, a c , 
is calculated separately in the RPA and is found to be 
a c — 6.4 a.u. The small valence-core (vc) a vc term that 
corrects the ionic core polarizability for the presence of 
the valence electrons is also calculated in the RPA; it is 
equal to —0.4 a.u. and —0.2 a.u. for the 6s 2 1 Sq and 
6s6p 3 Po states, respectively. DHF calculations are car- 
ried out as well for both of these contributions to evaluate 
the uncertainty associated with these terms, which was 
found to be negligible at the present level of accuracy. 

Accurate calculation of the polarizabilities of low- 
lying states is more difficult for Yb than for alkaline- 
earth atoms. It is known that the main contribu- 
tion to the ground state polarizability of Yb comes 
from 4f u 6s6p l P? and 4/ 13 5d6s 2 (7/2, b/2)° 1 states (see, 
e.g., [22|]). The energy difference between these states 
is only 3790 cm -1 and they strongly interact with each 
other. Calculations that treat Yb as an atom with only 
two valence electrons fail to account properly for the in- 
teraction between valence and core-excited states and de- 
scribe states with an unfilled / shell. While the state 
4/ 13 5d6s 2 (7/2, 5/2)° does not belong to the valence sub- 
space and is not directly mixed with the 4/ 14 6s6p l P° 
state in our calculations, its effect is introduced via the 
calculation of the effective Hamiltonian, since we allow all 
single and double excitations of the core shells during its 
construction. As a result, the polarizability calculation 
carried out via the solution of the inhomogeneous equa- 
tion does not appear to be affected by this problem. A 
theoretical explanation of this fact was suggested in [23[ 
which considered mixed and unmixed basis sets that 
included 4/ 14 6s6p 1 P 1 ° and 4/ 13 5d6s 2 (7/2, 5/2)° states. 
Excellent agreement of our results with all measured 
Yb polarizability-relatcd properties, including Stark shift 
and magic wavelength of the 6s 2 So — 6s6p 3 Pq transi- 
tion and the C§ coefficient of the Yb-Yb dimer, confirms 
that the mixing problem does not appear to affect such 
properties. We note that this is only true as long as no 
experimental data is substituted for theoretical quantities 
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TABLE I: The 6s 2 % and 6s6p 3 P ° static polarizabilities 
a g {uj = 0) of Yb and their difference Aa = a( 3 P °) - «( 1 5'o) 
calculated in CI, CI+MBPT, and Cl+all-order approxima- 
tions in a.u. The Cl+all-order values are taken as final. The 
present results are compared with other theoretical and ex- 
perimental values. 



Met no a 


aCSo) 


a( 3 PS) 


Aa 


CI 


187.9 


279.7 




CI+MBPT 


138.3 


305.9 




Cl+all-order 


140.9 


293.2 




Final 


141(3) 


293(10) 


152 


Theory [13] (2006) 


111.3(5) 


266(15) 


155 


Theory [22] (2007) 


143 






Theory [24] (2008) 


144.6 






Theory [23J (2010) 


141(6) 


302(14) 


161 


Ref. [25] a (2012) 


134.4 - 144.2 


280 - 290 




Expt. [11] (2012) 






145.726(3) 


"Constraints based 


on experimental data. The uncertainty 


in each of these values is 1.0. 





in any part of the calculations, since this will compromise 
the basis set completeness [23j . It follows that the direct 
solution of the inhomogeneous equation is expected to 
be more accurate than expected based on comparison of 
individual matrix elements with experiment. This con- 
clusion is important for future calculation of the parity- 
violating amplitudes that could by evaluated by the same 
techniques. 

Table U presents results for the static polarizabilities 
of the 6s 2 ^So and Qs6p 3 Pq states and their differences. 
We note that the states with an unfilled 4/ shell con- 
tributed less to the polarizability of the 3 Pq than 1 So 
state. In particular, even-parity states with an unfilled 
4/ shell lie rather high in energy and their contribu- 
tions to the polarizability and influence on other states 
is not so significant. The results obtained in the CI, 
CI+MBPT+RPA, and CI+all-order+RPA approxima- 
tions are presented. Our recommended values obtained 
at the CI+all-order+RPA stage are in a reasonable agree- 
ment with other theoretical values. We emphasize that 
our calculations are completely ab initio. The most re- 
cent recommended values of Ref. [23[ include adjustment 
to reproduce the experimental value of the magic wave- 
length. The set of accurate experimental data was used 
to set upper and lower bounds on the 1 So an d 3 Pq polariz- 
abilities in |25| . Our recommended values are in excellent 
agreement with these constraints taking an account the 
uncertainties. We can roughly estimate the uncertainty 
of our calculations as the difference of the CI+MBPT and 
Cl+all-order values, which yields 1.8% and 4.3% for 1 Sq 
and 3 Pq states. We note that the Cl+all-order value is 
higher than CI+MBPT for 1 S but lower for 3 P °, so we 
can expect that these uncertainties will add cumulatively 
for the Aa polarizability difference. However, our value 
of Aa agrees with a recent experiment to 4.3%, so our 
values are somewhat more accurate than the estimates 
above ( 1% and 3.5%, respectively). A direct measure- 



ment of the ground state polarizability with 1% accuracy 
would be an excellent test of the quality of calculations. 

To further check the accuracy of our approach we cal- 
culated the magic wavelength A for the ^Sq and 3 Pq states. 
At the magic wavelength, the frequency-dependent po- 
larizabilities of the two states are equal. We obtain 
A = 754 nm in the Cl+all-order approximation which 
is within 1% of the experimental value 759.355 nm [26j j . 
The polarizability of the 3 Pq state grows rapidly in the 
vicinity of the intersection of the ac polarizabilities. It 
means that even a small change in A leads to a signif- 
icant change in a( 3 P °). For example, the CI+MBPT 
value is significantly higher, 789 nm. Such close agree- 
ment of the Cl+all-order value with the experimental 
wavelength confirms the accuracy of the polarizabilities 
quoted above. 

An important application of the polarizability calcula- 
tion is to determine the shift of the Sq — 3 Pg transition 
frequency by the effects of the ambient blackbody radia- 
tion. The BBR shift is now one of the largest irreducible 
contributions to the budget of the uncertainty of opti- 
cal atomic clocks. The leading contribution to the BBR 
shift of the energy level g can be expressed in terms of 
its static polarizability a g (ui = 0) by [27j 



SE g = -^{a-Kf{k B Tfa g {Qi) [1 + r,] 



(2) 



where k B is the Boltzmann constant, T is the tempera- 
ture, and rj is a "dynamic" fractional correction to the 
total shift that reflects the averaging of the frequency de- 
pendence polarizability over the frequency of the black- 
body radiation spectrum. 

When the parameter \y n \ > 10, where Vn = {E n — 



Eg)/ '(k^T), rj can be be approximated by [27 



m + m + m 
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with 0.1% accuracy. We express the dominant term in 
the equation above as the second derivative of the polar- 
izability: 



m 



20 



21(2J 9 + 1) a fl (0) dE\ 



a g (0) 



(4) 



and find rji^So) = 0.00116 and ryi( 3 P ) = 0.00934. We 
calculated the second term in Eq. ([3]) using both a forth 
derivative of a and sum over states with the Cl+all-order 
values of the matrix elements; identical result r?2( 3 Po) = 
0.00029 was obtained. ry 2 is negligible for ^o, 0.000003. 
The third term can be neglected at the present level of 
accuracy. The resulting values of the dynamic corrections 
to the BBR shift at 300 K are Svbbr^Sq) = -0.0014 Hz 
and <5^bbr( 3 Po) = —0.0243 Hz, respectively. 

The total dynamic correction to the BBR shift at 
300K is determined as the difference between the indi- 
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TABLE II: The values of the Ce coefficients (in a.u.) for 
the homonuclear Yb dimer and the heteronuclear alkali- 
metal/group II - Yb dimers. All atoms are in their ground 
states. The a(iw) for alkali and Mg, Sr, and Ca are taken from 
Ref. [H in rows CI+MBPT and CI+all (a) . The a(iuj) for Mg, 
Sr, and Ca are calculated with the Cl+all-order method in the 
present work in row CI+all (b) . The present CI+MBPT and 
Cl+all-order Yb a(iu) values are used in rows CI+MBPT 
and Cl+all-order'"' 6 '' , respectively. The uncertainty of the 
final Cl+all-order values is estimated to be 2%. 





Li-Yb 


Na-Yb 


K-Yb 


Rb-Yb 


Cs-Yb 


CI+MBPT 


1534 


1655 


2548 


2807 


3367 


CI+all (a) 


1551 


1672 


2576 


2837 


3403 




Yb-Yb 


Mg-Yb 


Ca-Yb 


Sr-Yb 




CI+MBPT 


1901 


1086 


2000 


2414 




CI+all (a) 




1093 


2017 


2435 




CI+all (b) 


1929 


2024 








Expt. [29] 


1932(35) 











vidual shifts, <5f BB p = —0.0229 Hz. Combining this re- 
sult with the experimental determination of the Acv = 
145.726(3) a.u. [Il|, we get the final result for the 
BBR shift at 300 K: £z/ BBR = -1.2777(8) Hz. This 
value is in excellent agrement with the determination 
of the BBR shift mostly from the experimental data 
^bbr = —1.2774(6) Hz that was just reported in [30| . 
We have verified that the Ml and E2 contributions to the 
BBR shift are negligible at the present level of accuracy. 
Details of the calculation of the dynamic correction to the 
BBR shift are given in the supplementary material (20| . 

Many of the same considerations concerning accu- 
rate calculation of the frequency-dependent polarizabil- 
ity arise in the calculation of the van der Waals coeffi- 
cients. If two atoms A and B have spherically symmet- 
rical ground states, the leading power of the long-range 
interactions takes the form V(R) = —Cq B /R 6 , where R 
is the distance between atomic nuclei. The the van der 
Waals coefficient Cq B can be calculated as [HI 

C^ B = - a A {iw)a B {iu)duj, (5) 
n Jo 

where a(iu>) is the frequency-dependent polarizability at 
an imaginary frequency. In practice, we compute the 
Cq B coefficients by approximating the integral (O by 
Gaussian quadrature of the integrand computed on the 
finite grid of discrete imaginary frequencies |32| . 

For the alkali and group II atoms, we use frequen- 
cies and weights tabulated in Ref. !28j at 50 points. 
These dynamic polarizabilities were obtained by com- 
bining high-precision experimental data for matrix ele- 
ments of principal transitions with high-precision many- 
body methods, such as linearized coupled-cluster ap- 
proach and CI+MBPT. The accuracy of the correspond- 
ing homonucelar Cq was estimated to be better than 1% 
for all cases relevant in this work with the exception of 
Ca, were it was 1.5% (281 ] . 



The Yb imaginary frequency polarizabilities Q/iu) for 
the ground state are calculated in this work by solv- 
ing the inhomogeneous equation (JlJ with the appropri- 
ate modifications. We use the same 50-point frequency 
grid as in Ref. [28[ for consistency. To evaluate the 
uncertainty in the Cg coefficients, we carried out both 
CI+MBPT and Cl+all-order calculations of the ground 
state a(iu>). The same alkali and group II data are used 
in both cases. The results are summarized in Tablellll We 
find that the differences between CI+MBPT and Cl+all- 
order results are actually smaller (1-1.5%) than for the 
ground state static polarizability (1.8%) since the differ- 
ences decrease with u> for a(ioj). As a result, we expect 
the accuracy of the Cq to be on the same order as the 
static polarizability, rather than larger by a factor of two. 
Moreover, our value for the Cg coefficient of the homonu- 
clear Yb dimer is in excellent agreement with the exper- 
imental result [29(, which is accurate to 1.8%. Based on 
the comparison of the CI+MBPT and Cl+all-order val- 
ues for heteronuclear Cq coefficients, and agreement with 
experiment for the Yb 3 Pq — 1 Sq Stark shift and magic 
wavelengths, and Ce coefficient for Yb-Yb dimer, we es- 
timate that our predictions of the C§ coefficients for the 
heteronuclear alkali- metal atom/group II - Yb dimers are 
accurate to about 2%. 

In conclusion, we have carried out fully ab initio all- 
order calculations of Yb properties. Our values of the Yb 
— 3 Pq Stark shift and magic wavelength as well as the 
Cq coefficient of the Yb2 dimer are in excellent agreement 
with experiment. We have developed a new approach of 
calculation of the dynamic correction to the BBR shift 
that does not involve an explicit sum over states. The 
Stark shift of the clock transition was determined exper- 
imentally [ll[ with a high precision. As a result, the un- 
certainty in the dynamic correction can now be directly 
related to the uncertainty of the BBR shift of this transi- 
tion. Thus, when combined with the recent measurement 
of the Yb clock Stark shift [ll|, our calculation of the 
dynamic correction allows us to reduce the fractional un- 
certainty due to the BBR shift in the Yb optical lattice 
clock to 10~ 18 level. The same method can be used to 
evaluate the dynamic correction for any optical atomic 
clock. Finally, we have presented the first recommended 
values of Ce coefficients for alkali/group II- Yb dimers for 
future experimental efforts in producing ultracold polar 
molecules. 
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